{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explain LSTM->XGB on TBI signal data\n",
    "\n",
    "Explain XGBoost model trained on LSTM embedded data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"4\"\n",
    "\n",
    "import numpy as np\n",
    "from tbi_downstream_prediction import split_data\n",
    "\n",
    "import keras\n",
    "from keras.utils import multi_gpu_model\n",
    "from keras.layers import Input, LSTM, Dense, Dropout\n",
    "from keras.models import Sequential, load_model, Model\n",
    "from matplotlib import cm, pyplot as plt\n",
    "from sklearn import metrics\n",
    "from os.path import expanduser as eu\n",
    "from os.path import isfile, join\n",
    "import numpy as np\n",
    "import random, time\n",
    "\n",
    "import tensorflow as tf\n",
    "from keras.backend.tensorflow_backend import set_session\n",
    "config = tf.ConfigProto(allow_soft_placement=True,gpu_options = tf.GPUOptions(allow_growth=True))\n",
    "set_session(tf.Session(config=config))\n",
    "\n",
    "PATH = \"/homes/gws/hughchen/phase/tbi_subset/\"\n",
    "DPATH = PATH+\"tbi/processed_data/hypoxemia/\"\n",
    "data_type = \"raw[top11]\"\n",
    "\n",
    "feat_lst = [\"ECGRATE\", \"ETCO2\", \"ETSEV\", \"ETSEVO\", \"FIO2\", \"NIBPD\", \"NIBPM\", \n",
    "            \"NIBPS\",\"PEAK\", \"PEEP\", \"PIP\", \"RESPRATE\", \"SAO2\", \"TEMP1\", \"TV\"]\n",
    "\n",
    "# Exclude these features\n",
    "weird_feat_lst = [\"ETSEV\", \"PIP\", \"PEEP\", \"TV\"]\n",
    "feat_inds = np.array([feat_lst.index(feat) for feat in feat_lst if feat not in weird_feat_lst])\n",
    "feat_lst2 = [feat for feat in feat_lst if feat not in weird_feat_lst]\n",
    "\n",
    "y_tbi = np.load(DPATH+\"tbiy.npy\")\n",
    "X_tbi = np.load(DPATH+\"X_tbi_imp_standard.npy\")\n",
    "\n",
    "X_tbi2 = X_tbi[:,feat_inds,:]\n",
    "(X_test, y_test, X_valid, y_valid, X_train, y_train) = split_data(DPATH,X_tbi2,y_tbi,flatten=False)\n",
    "\n",
    "PATH = \"/homes/gws/hughchen/phase/tbi_subset/\"\n",
    "RESULTPATH = PATH+\"results/\"\n",
    "label_type = \"desat_bool92_5_nodesat\"\n",
    "lstm_type = \"biglstmdropoutv3_{}\".format(label_type)\n",
    "RESDIR = '{}{}/'.format(RESULTPATH, lstm_type)\n",
    "if not os.path.exists(RESDIR): os.makedirs(RESDIR)\n",
    "GPUNUM = len(os.environ[\"CUDA_VISIBLE_DEVICES\"].split(\",\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load model and generate independent explanations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import shap\n",
    "import xgboost as xgb\n",
    "\n",
    "from tbi_downstream_prediction import *\n",
    "PATH = \"/homes/gws/hughchen/phase/tbi_subset/\"\n",
    "DPATH = PATH+\"tbi/processed_data/hypoxemia/\"\n",
    "RESULTPATH = PATH+\"results/\"\n",
    "\n",
    "# Load data\n",
    "X_train = np.load(DPATH+\"X_train_embed_lstmbigdropoutv3_50n_arr.npy\",mmap_mode=\"r\")\n",
    "X_valid = np.load(DPATH+\"X_valid_embed_lstmbigdropoutv3_50n_arr.npy\",mmap_mode=\"r\")\n",
    "X_test  = np.load(DPATH+\"X_test_embed_lstmbigdropoutv3_50n_arr.npy\",mmap_mode=\"r\")\n",
    "\n",
    "# Set important variables\n",
    "label_type = \"desat_bool92_5_nodesat\"; eta = 0.02\n",
    "hosp_data = \"tbi\"; data_type = \"lstm_big_50n[top11]\"\n",
    "mod_type = \"xgb_{}_eta{}\".format(label_type,eta)\n",
    "\n",
    "# Set up result directory\n",
    "RESDIR = '{}results/{}/'.format(PATH, mod_type)\n",
    "if not os.path.exists(RESDIR): os.makedirs(RESDIR)\n",
    "\n",
    "# Set parameters to train model\n",
    "param = {'max_depth':6, 'eta':eta, 'subsample':0.5, 'gamma':1.0, \n",
    "         'min_child_weight':10, 'base_score':y_train.mean(), \n",
    "         'objective':'binary:logistic', 'eval_metric':[\"logloss\"]}\n",
    "\n",
    "save_path = RESDIR+\"hosp{}_data/{}/\".format(hosp_data,data_type)\n",
    "bst = xgb.Booster()\n",
    "bst.load_model(save_path+'mod_eta{}.model'.format(param['eta']))\n",
    "X_tbi2 = X_tbi[:,feat_inds,:]\n",
    "(_, _, _, _, X_train_raw, _) = split_data(DPATH,X_tbi2,y_tbi,flatten=False)\n",
    "\n",
    "# Generate explanations\n",
    "import numpy as np\n",
    "np.random.seed(10)\n",
    "background_inds = np.random.choice(np.arange(0,X_train.shape[0]),1000)\n",
    "X_train = np.load(DPATH+\"X_train_embed_lstmbigdropoutv2_arr.npy\")\n",
    "X_train_background = X_train[background_inds]\n",
    "ind_explainer = shap.TreeExplainer(bst,data=X_train_background,feature_dependence=\"independent\")\n",
    "X_train_ind_explanations = ind_explainer.shap_values(X_train)\n",
    "np.save(save_path+\"X_train_explanations_ind\",X_train_ind_explanations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the aggregated explanations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_ind_explanations_agg = np.transpose(np.vstack([X_train_ind_explanations[:,i*50:(i+1)*50].sum(1) for i in range(len(feat_lst2))]))\n",
    "shap.summary_plot(X_train_ind_explanations_agg,X_train_raw.sum(1),feature_names=feat_lst2)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
